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1 Introduction 



Mass-loss due to various processes during the lifetime of a star will increase the semi-major axis of a 



wide binary (Dommanget, 1963; Valls-Gabaud, 1988), what we refer to here as orbital separation 
amplification. It is expected that mass-loss during the stellar lifetime, in particular Post-Main 
Sequence (PMS) mass-loss, will statistically distort a frequency distribution of weakly interacting 
wide binary (Fragile Binary, FB) separations. Visual confirmation of such processes would provide 
tests of current theories of stellar evolution and fundamental information on the dynamics of 



Galactic structures (Wasserman, 1988). 



As discussed by Oswalt et al. (2007) and Chaname (2007), FB's provide insight into several 
astrophysically important problems that cannot be easily addressed by other means. In particular, 
such fragile pairs provide valuable constraints on the Galactic tidal potential, stellar and ISM 
perturbations, and PMS mass-loss, all of which will affect orbital separations. A number of early 



papers have addressed this topic: Eggen ( 1965 ) studied the colors and magnitudes of 228 visual 



binaries (VBs), deriving an vs. R — I calibration from trigonometric parallaxes. Van Dessel 



and Dommanget (1980) showed that, using proper motions from micrometer measurements, a pair 



could be deemed gravitationally bound and thus, through statistical studies and interferometry. 



orbital parameters could be determined. Trimble (1974) performed a statistical study of the 
observed populations of binary systems (798 VBs) and was one of the first to note the peak in the 
mass ratio (q) distribution at g ~ 1. 

In recent years, computing power has improved to the point that modeling the evolution of 



whole stellar binary populations is now possible. Examples include: Wasserman (1988); Valls- 



Gabaud (1988), and more recently Hurley et al. (2001). Prior work has shown evidence for orbital 



separation amplification by as much as half an order of magnitude in binaries, where at least one 



of the partners has undergone PMS mass loss (Hadjidemetriou 1966 Greenstein 1986 Wood and 



Oswalt, 1992). Considering that the majority of all stars are born with one or more companions 



as a result of the initial conditions required for a stable system (Mathieu, 1994), the importance 
of understanding how mass-loss affects the evolution of binary systems cannot be over stated. 
In addition, the initial-to-final mass relation (IFMR) for WD stars is at present very poorly 



constrained (Weidemann, 1993; Zhao et al. , 2011B). Proper modeling of FB orbital separation 
amplification can provide a valuable independent estimate of this important relation. 
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The distribution of WD + MS pair separations is most strongly affected by the PMS (i.e. red 
giant phase) mass-loss. FB systems typically have component separations of ~ 10^ AU (Zinnecker 



1984); this results in negligible tidal interactions and mass transfer between companions. Such 



pairs evolve as if they are two separate but coeval stars. Once mass-loss has ended, orbital 
evolution is limited to the effects of external forces such as stellar encounters, giant molecular 



cloud encounters, gravitational tidal forces, and rare encounters with other stars (Chaname and 



Gould 2004 Jiang and Tremaine, 2010) . 



In this study, which updates and improves upon our previous work (Johnston et al. , 2007), we 
compare a theoretically generated population of FB angular separations to the observed angular 
separation distributions of FBs found in a large-scale survey. We constructed a simulation code, 
which we refer to here as the Double Star (DS) program, to generate the initial parameters and 
distributions of the tested population. The DS program evolves the binary components using the 



Single Star Evolution (SSE) code (Hurley et al. , 2000A) over the lifetime of the simulation. The 
DS program is also responsible for the evolution of individual orbital elements and translation 
of the modeled physical properties into observable features. The results of the translation is a 
population of theoretical FBs, with observable properties. 



2 Observational Data 



It should be possible to produce a statistically significant frequency distribution of angular sepa- 
rations for FBs with various constituents, by sampling a large survey such as the Luyten Double 



Star (LDS) Catalog (Luyten, 1969). While the LDS catalog, with 6124 pairs and multiples, does 
not have the volume of more recent surveys such as the Sloan Digital Sky Survey (SDSS), it 
provides a well-observed large sample to compare to our theoretical models. The identification 
of FB candidates requires selection of an upper bound separation. Luyten found pairs separated 
by as much as 400". Due to the POSS survey plate scale, Luyten's minimum angular separation 
limit was ~ 0.5". We will use this sample as a first test of our programme, while we note that 



further samples of fragile binaries are currently being constructed (e.g., Sesar et al. 2008 Quinn 



and Smith 2009 Longhitano and Binggeh 2010; Zhao et al. 2011A) 
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3 Initial Parameters 



To build a synthetic model of the local Galactic FB population, we first proposed a set of initial 
parameters that defined the individual components of our binaries, as well as the binary system 
separations. We used a combination of published functions from previous theoretical models and 
observational surveys. When there were no such references, or when the functions published were 
not appropriate for FBs, we implemented physically reasonable distributions. 



3.1 Mass Generating Function 



Following Kroupa et al. (1993) we assumed a Salpeter distribution of ZAMS stars using a variation 



of Eggleton's mass generating function (Equation [T|. 



M{x) = 0.33( 



1.0 



1.0 



:i.o 



XI 



(1) 



;i.0-x)0-75 + 0.04(1.0 1.04 
Where x is a uniform randomly distributed number and M(x) is the randomly generated, 
restricted, mass of a ZAMS star. For this project, we initially assumed a minimum mass of 



O.OTMq (spectral type-T), given by Chabrier (2003). From this we produced an initial distribution 
of primary masses (Figure [T]). 



3.2 Mass Ratio 



The mass ratio, q, is one of the most poorly determined variables in a wide binary system. Current 



observed distributions fall into two categories: (1) a distribution that peaks at g ~ 1 (Trimble 



1987) or (2) a bimodal distribution with peaks at g ~ 0.2 and q ~ 0.7 (Halbwachs et al. , 2003). 
The former are explained in terms of turbulence and fragmentation during the binary formation 



process (Lucy and Ricco, 1979). This tends to result in near equal clumps of would-be stars forming 



out of a single GMC Bimodal distributions seem to be created by selection effects (spectroscopic 
methods vs. common proper motion) and evolutionary effects (binary capture, disruption, mass- 
loss, mass-transfer,etc.), for references of either see Kruszewski ( 1967 ) ; [Paczynski (1971); Shu and 



Lubow (1981). For the purpose of finding an initial mass ratio function (IMRF), we will tend to 



use the discussion of distributions which tend towards a maximum of one in the set 0< q <1 (unity- 
like), as these will be the most likely candidates to explain a mass-ratio distribution created by 
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Figure 1: Histogram of the distribution of primary masses from initial generation (in red with 45 
deg. hatching) and after evolution and orbital expansion (in blue with -45 deg. hatching). The 
difference in population size is resulting from the dissolution of some binaries from the original 
population. Each bin is of size 0.05 in log (Mq) 

a GMC to stellar binary creation process. After multiple trials with various models, we chose a 
random uniform distribution, where 0< q <l and A is a constant (Equation [2]). 



fiq) = A (2) 

A constant mass ratio distribution function has the added advantage of producing an undis- 
torted secondary mass distribution, implying that the same physical processes produce the primary 
and secondary components. Using an acceptance-rejection method of random variate generation 



(Law et al. , 2005), we produced the distribution of secondary masses seen in Figure |2j In addi- 
tion we restricted all lower mass components to Mgtar > O.OTM©. This has the added effect of 
producing a tendency for g ~ 1, as seen in Figure |3j in accord with prior observational studies. 
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Figure 2: Histogram of the distribution of secondary masses from initial generation (in red with 
45 deg. hatching) and after evolution and orbital expansion (in blue with -45 deg. hatching). 
The difference in population size resulted from the dissolution of some binaries from the original 
population. Each bin is of size 0.05 in log (Mq) 

3.3 Orbital Parameter Distribution Functions 

In this paper, we focused on the development of two specific orbital element distributions, specif- 
ically the distribution of semi-major axis (separation) and the distribution of eccentricity. From 
these parameters, we develop elliptical binary mechanics for each binary in our initial population. 
The statistical characteristics of the separation distribution for the un-evolved binaries, as well as 
the distribution characteristics for the primary masses, secondary masses, and periods is given in 
Table [l] at the end of this section. 
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Figure 3: Histogram of the distribution of mass ratios from initial generation (in red with 45 deg. 
hatching) and after evolution and orbital expansion (in blue with -45 deg. hatching). The differ- 
ence in population size resulted from the dissolution of some binaries from the original population. 
Bin size is of size 0.05 in g-ratio 

3.3.1 Initial Separation Function 

The separations of binary systems is a difficult observable parameter, requiring distance, incli- 
nation, eccentricity and mass. One of the first and most often used references to separation 



distribution was given by Opik (1924), and implemented by Zinnecker (1984) and Poveda (1988). 



It is a simple uniform log distribution: 



/(a) = log(a) 



(3) 



Where a > 100 A.U., approximately the closest separation where post-MS mass exchanges are 
unlikely to have occured. We transform equation |3] into a random variate separation generating 
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function (Equation |4]). 

a(X) = 1020^+20 (4) 

Where < X < 1 and X is a uniformly distributed random number. Similarly, we used the 
separation and masses now assigned to each pair, along with Kepler's 3'^'^ law to calculate the 
period for each binary. 




Figure 4: Histogram of the distribution of separation from initial generation (in red with 45 deg. 
hatching) and after evolution and orbital expansion (in blue with -45 deg. hatching). The differ- 
ence in population size resulted from the dissolution of some binaries from the original population. 
Each bin is of size 0.05 in log {A.U.) 

In addition to the physical limit of a > lOOA.U. for non-interacting pairs, we also imposed 



observed limitations on our binary pairs. We rely on the Zinnecker (1984) definition of CPMB/FB, 
VB, and SBs and limit the initial smallest separation that a common proper motion binary can have 
as lOOA.U.. This is wide enough that the two components of the pair are individually resolvable 
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Figure 5: Histogram of the distribution of periods from initial generation (in red with 45 deg. 
hatching) and after evolution and orbital expansion (in blue with -45 deg. hatching). The differ- 
ence in population size resulted from the dissolution of some binaries from the original population. 
Each bin is of size 0.1 in log [Myr) 

and the effects of Roche Lobe interaction and component-to-component tidal interactions are 
minimal compared to Galactic effects and stellar GMC encounter perturbations. 

We additionally imposed a maximum of 10^ A.U. , which appears to be the maximum theoretical 



separation in wide pairs (Wasserman and Weinberg, 1987) imposed by gravitational disruptions 



of the Galaxy and stellar perturbations. An algorithm was developed to impose this upper limit 
continuously over the evolution time frame and will be discussed below in section 4.3. 

3.3.2 Eccentricity 



The simulation derives the initial distribution of eccentricities following Heggie (1975). The initial 
distribution of eccentricity for wide binary pairs is given in Equation |5] 
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/(e) = 2e (5) 
We transform this distribution into the random variate generating function (Equation [6]). 



e(X) = (6) 

Where < e < 1 and XisO<X<l,a uniformly distributed random number. This 
distribution was derived based on assumptions of a Maxwelhan particle distribution. It becomes 
apparent from the generating function that binaries, with this eccentricity distribution, are not 
favored to have initial circular orbits. In fact, using this distribution one can easily see that the 
expected initial density of circular orbits is zero. For wide binary pairs we can convince ourselves 
that this is true, as these pairs tend to have coeval components but have had little or no interaction 
beyond the gravitational binding which makes them a pair (as discussed in the mass ratio section). 
This creation via loose interaction would tend to result in more initially eccentric orbits. 

Table 1: Statistical Characteristics of the Un-evolved Distributions 



Metric 


Average 


StdDev 


Skew 


Kurtosis 


Min 


Max 


Mass 1 (log Mo) 


-0.28 


0.41 


0.58 


0.79 


-1.10 


1.93 


Mass 2(logM0) 


-0.55 


0.40 


0.92 


1.03 


-1.10 


1.84 


Separation {log A.U.) 


3.00 


0.58 


-0.002 


-1.20 


2.00 


4.00 


Period (log Myr) 


-1.46 


0.826 


-0.008 


-1.08 


-4.04 


0.39 



3.4 Birth Rate 



Our experiments incorporated an initial stellar burst of 10^ binaries with Willems and Kolb (2004) 
constant repeating burst model of the stellar birth rate, 10^ binaries being created every lOMyrs, 
resulting in a total initial population of 100,100 binaries. All bursts have identical initial factors 



and distributions, resulting in a continuous and smooth distribution (Hurley et al. , 2000A) in the 



initial population. The only references that challenge this approach, suggest that metallicity may 
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play a role in distorting the initial distribution function. However, the region of space that we 
are attempting to model is limited to the solar neighborhood (i.e. distances under lOOpc), FBs in 
this region tend to have solar metallicitys. We therefore assumed solar metallicity for all ZAMS 
stars. Future work may address a more realistic history of star formation rate, dervied for the 



solar neighbourhood (Hernandez et al. , 2000). 



3.5 Color Estimation 

The goal of this project is to compute a theoretical population of FBs that accurately compared to 
observed angular separation distributions. This requires a transformation of our physical charac- 
teristics for each binary system into observables such as color, magnitude and angular separation. 



We employed a mix of Kawka and Vennes (2006); Kurucz (1979); Hawley et al. (2002) temperature 



vs. color relations for MS and WD stars, to produce the ugriz{J) absolute magnitudes required 
for comparison. 

We divided stars into four categories: MLT, WD, MS, and other. MLT refers to the spectral 
type of the star: M, L or T, corresponding to effective temperatures from ~ 3800/^ to lOOO-ft" 



(Reid, 2000, 2002). The spectral type of a star below 3400/^ is assumed to be dependent only on 



temperature. We justify this estimation with three points: (1) Reliable statistics are not available 
for surface gravities of low mass stars (early-T and late-L); (2) The current classification scheme 
for MLT dwarfs is based mostly on line features. A number of flux ratios {H2O, CO2, CH^^^etc) 



show linear relations with spectral type (Burgasser, 2002); (3) Stephens (2001) used a hybrid 
approach to show that temperatures and spectral types scale linearly with each other for L-type 
dwarfs. 

Next, we used a combination of spectral type and temperature relations to transform the 



observed absolute magnitudes by Hawley et al. (2002) in r*i*z*J to produce a temperature vs. 



absolute magnitude relation. We then used this relation to produce colors and errors for all MLT 
stars. Neither the LDS nor the SDSS survey provides apparent J magnitudes. Future surveys 
however will produce near-infrared colors (y). At such a time our code can be updated and 
expanded. 



For WD and MS stars, we use models provided by Kawka and Vennes (2006) and Kurucz 



(1979), respectively. These atmospheric models are grids of color as functions of temperature 
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and gravity. We employed a 2nd/3rd order unevenly spaced 2-dim interpolation subroutine to 
determine both expected absolute magnitude and error in absolute magnitude given T and \og{g). 
As before, we estimated u*g*r*i*z* colors. This allowed simple estimates of distance modulus for 
each pair. 

Models used to compute discrete grids are availably derived from the Kurucz stellar atmo- 
sphere models. For MS stars temperatures range from 3, 170K <T< 42, OOOK and < log{g) < 
5. For WD stars, two sets of grids were used: a fine resolution with range 7, OOOK < T < 84, OOOi^ 
and 7 < log{g) < 9.5 in 0.5 intervals and a broad resolution set with 4, 5007^^ < T < 84, OOOK and 
7 < log(5') < 9 with 1.0 gravity intervals. The fine resolution set was used for stellar temperatures 
with 7, OOOK < T < 84, OOOK. For temperatures above 84, OOOK and below 7, OOOK, the broad 
resolution set was used. 



3.6 Distance 

In order to compute the distance modulus, a variety of models were tried, including a synthetic 



Galactic population model that includes distance as a function of mass (Robin et al. , 2003). 
Initially, we produced a distance model that best fits the observed LDS distribution with the 
given parameters. The adopted distance distribution had the form: 

f{d) = log(rf) (7) 
This equation was then transformed into a random variate separation generating function: 



a{X) = 10 



3.0X+1.0 



(8) 



Where < X < 1 and X is a uniformly distributed random number. We expect this will, 
through random evolutionary processes, produce the synthetic lognormal distribution of distances 



that Robin et al. (2003 ) obtained. We also assume that binary and single star distance distributions 
are identical. 



^http://www. stsci.edu/hst/observatory/cdbs/k93models. html 
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3.7 Stellar Evolution Code 



One of the two major tasks in this project is the abihty to accurately evolve a theoretical star while 
keeping track of its mass, luminosity, and temperature; as previous analysis used analytical, rather 
than physically-motivated, mass- loss rates (e.g. Dommanget||1963 Docobo et al.||1998 Rahoma et 



al. 2009). For a population of binary stars this requires substantial computational resources. As 



mentioned above, the Single Star Evolution (SSE) program (Hurley et al. , 2000A) was adapted 
for this effort. The basis of the SSE is a set of comprehensive analytical formulae, that accurately 
track the evolution of a star with given initial mass and metallicity. The goal of the SSE program 
was to allow for the rapid calculation of a number of stellar properties for a given stellar lifetime, 
thus facilitating the production of a population synthesis model that could produce a statistically 
significant sample of binaries in a reasonable computational timeframe. 

The key to the SSE code's accuracy and speed is its process style, which employs a combina- 
tion of interpolation formulae and tabular construction. The HPT group uses tabulated stellar 



1998 


)• 


Pols et al. 


(1998 



later expanded the range of mass and metallicity of the tables incorporated in the SSE program. 
Analytic approximations to the movement of a star in the Hertzsprung-Russell diagram (HRD) 



are a simplistic approach to a complex problem (Eggleton et al. , 1984). However, supporting these 



calculations with tabulated stellar models allows a more precise estimation of physical properties 
and greatly decreases computational time. 

The SSE does not take into account stars with solar masses lower then O.2M0. For stars with 



masses lower then this limit we adopted the Baraffe et al. (1998) low-mass stellar evolutionary 



models. We do not expect such low mass stars to undergo stellar mass evolution, but we do expect 
small changes in other stellar parameters over time. Using these evolutionary models we were able 
to determine the initial and final values of temperature and gravity based on initial mass and 
metallicity. 

We selected the following parameters for the SSE code: no wind-loss, a Reimers Coefficient for 
mass-loss of 0.5, and a constant metallicity distribution of 0.1. This set of parameters provides an 
initialization point from which to develop further analysis, while still generating realistic results. 
Future research may include analysis (ANOVA) to determine of optimal parameter distributions 
and values for comparison of theoretical and observed population distributions. The authors of 
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the SSE code have subsequently developed a follow-up algorithm, the Binary Star Evolutionary 
code. For development purposes, this new program was not used, as the primary concern of the 
DS-code is the generation of evolved binary populations and information, relating to fragile coeval 
binaries. The authors feel that wrapping the elliptical orbital computations, stellar population 
generating functions, and apparent magnitude translations around the original SSE code, allows 
for the maximum amount of computational efficiency and provides the necessary forensic test 
points for further experimentation, specifically in the generation of an ANOVA experiment for 
optimal parameter and distribution selection. 



3.8 Orbital Evolution 



In the present model, we assumed isotropic non-conservative mass-loss for both stars, both during 
the MS and PMS lifetimes. Since we modeled what could be considered two phases of mass-loss, 
slow (MS) and sudden (PMS), we found a two stage evolutionary code was needed. Following 



Valls-Gabaud (1988); Hadjidemetriou (1966); Kopal (1978), we considered the effect double, non- 



interacting, mass-loss, has on the orbital parameters of the binary. For this experiment we modified 



versions of Valls-Gabaud (1988) equations into the simulation code presented here. We allowed 



no accretion onto either star from its companion and no accretion from the ISM or any outside 
bodies. The addition of perturbation forces from exterior sources beyond those addressed within 
this work, such as loose tertiary components, dense stellar populations, or other chance encounters 



were not modeled here, however we would recommend Heggie (1996) as starting research in the 
matter. 

We can show via adjustment of the mass- loss orbital parameter equations that the perturbing 
force due to an isotropic variation of mass is tangential to the orbit, thus only the planar orbital 
elements will be perturbed. This assumption holds when the mass-loss timescale of the binary 
is much greater than the orbital timescale (period). These perturbation equations are given as a 
function of time, the logarithm of the total system mass, the eccentric anomaly, and the eccentricity 



are given as the equations M and 10 



d log a 1 + e cos E d log M 



dt 



1 — e cos E dt 



(9) 
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de 
di 



e^) * cosE dlogM 



1-ecosE dt ^^^"^ 

Where is the rate of mass- loss over the time period of interest (iteration), and E is the 

eccentric anomaly. Note, that when expressing these equations in terms of the eccentric anomaly, 

Kepler's equations is no longer valid and should be replaced with the equation. 



dE 
~dt 



27r dlogMsinE 

\ 

P dt e 



[1 — e* cosE) ^ 
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The effect that this MS mass-loss has on the binary system is documented by Verhulst (1969), 



Dommanget (1963), Docobo et al. (1998) and Rahoma et al. (2009). Here, Verhulst focused the 



majority of the research on the effect a Jeans-Eddington mass-loss relation would analytically have 
on the eccentricity of the binary. It was shown that one can expect periodic variation about a 
mean eccentricity for a binary with isotropic continuous mass-loss depending on the parameters of 



the Jeans-Eddington mass-loss relationship. Hadjidemetriou (1963) demonstrated similar effects 



showing that with decreasing system mass the eccentricity of the system will sinusoidally increase. 



More recently Veras et al. (2011 ) demonstrated orbital parameter evolution for low q-ratio binaries 



(g < 0.1, exoplanets and smaller bodies). They showed that, at least for these conditions, the 
effect of mass-loss, varies greatly depending on the initial eccentricity distribution. 

Given that our mass-loss is a function of the stellar evolution code designed into the SSE, 



with the additional consideration of the low mass evolution models provided by (Baraffe et al. 



1998), we expect a similar result (increasing non-linear, sinusoidal eccentricity as a function of 



mass-loss). However we expect that the numerical (and probabilistic) nature of our experiment 
would demonstrate a different (and potentially more complete) relationship. For situations where 
mass-loss occurs on a timescale less than that of the period of the binary, we follow the sudden 



mass-loss equations 12 and 13 outlined by Hills (1983). 




AM 

M 



'ao\ AM 
. r ) M 



1-(1 



el)* 



2ao AA/ 
r M 



n \ 0.5 



^1 _ AM\2 



(12) 



(13) 



We have further transformed the original equations into a function of the mass-loss ratio over 
a given defined time period (delta mass over original mass). The increase in separation and 
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eccentricity, is a function of how close to periapsis or apoapsis the binary is at the time of mass- 
loss. It can be shown that for binary separation distances less than the semi-major axis (which is 
dependent on the eccentricity); sudden mass-loss can drive the binary together (decrease orbital 



separation). Similarly, as stated by Hills, if mass-loss occurs in a binary with a very eccentric orbit, 
the average post-explosion eccentricity will be less than the pre-explosion eccentricity. Conversely, 
given the position of the binary, the separation can be greatly amplified, resulting in populations 
of WD + MS and WD + WD pairs with separations on average larger than their MS + MS 
counterparts. We can demonstrate for given eccentricity and location in the orbit, post-main 
sequence mass-loss can result in a perturbation of the orbital eccentricity beyond a bound elliptical 



orbit, i.e., e > 1. This is demonstrated in Equation 14 



r AM , , 

Ya < -Jf 

It is apparent, that increasing eccentricity can pose a problem for the binary system. The 
potential for binary separation can result from increasing the eccentricity of the binary beyond 
(e > 1). Thus, we expect that the increase in eccentricity, and therefore the associated orbital 
velocities, in combination with the wide orbital amplitudes, will result in a number of binary 
systems being perturbed out of their orbits (dissolution via external perturbation). It should be 
noted that decreases in eccentricity resulting from sudden PMS mass-loss, coupled with the lack 
of mass-loss associated with white dwarf components, can result in the effective circularization of 
the orbit. This is apparent from the disproportionate number of binaries in the WD + MS and 
WD + WD population evolving to a final eccentricity under 0.1 in our simulations. 

To model the effect high eccentricity has on the binary system, specifically on the dissolution 
rate of the binary populations resulting from isotropic mass-loss and galactic tidal forces, we have 
implemented a function that fiags evolved binary systems which have eccentricity greater than 



0.95. This, combined with the Antonov (1972) tidal limit provides a first order approximation 



to dissolution effects on the binary systems. Other improvements are outlined in the Binary 
Dissolution section 
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3.9 Binary Dissolution 

Even a simple model for the orbital evolution of wide binaries must include the influences of stellar 



encounters, Giant Molecular Clouds (GMCs) and the Galactic tidal field (Weinberg, 1988). For 



the wide separations of FBs, even minor forces such as these operating over the course of millions 



of years, can affect the binary separations, eventually dissolving them. Antonov (1972) described 



a tidal limit for wide binaries, where the gravitational forces from Galactic tidal perturbations 



outweigh the gravitational attraction between the binaries themselves (Equation 15). 



ar ~ 2.07 x IQ^A.U. 



Mrs 



(15) 



This relation does not include either GMC or stellar encounters, but it does provide an initial 



cut off for galactic tidal disruption as a function of mass. We take the Antonov (1972) equation 



as a first approximation to a dissolution algorithm. A more realistic approach would include 
all disruption effects. A rigorous treatment is beyond the scope of the current project. For a 



first-order approximation to the dissolution time problem, we used the Wasserman and Weinberg 



(1987) Galactic tidal perturbation formula for a 2Mq binary system with q = 1. The expected 



half- life (^1/2) of an FB system with no consideration for the GMC encounter is given by Equation 



ti/2 = 4.3 X lO^Myr 



-1.4 



(16) 



^2.07 X 10^^ 

For wide binaries with consideration for the GMC local density of 3.6 x lO^pc"^, this relation 



becomes Equation 17 



ti/2 = 1.32 X 10^ Myr 



^1.4 



(17) 



^2.07 X 10^^ 

We define ti/2 as the 50% likelihood survival time for a binary with separation a in A.U. 
While this cannot be used to estimate the expansion caused by a GMC and perturbations, using 
a binomial distribution selection program we produced a dissolution algorithm to remove binaries 
from the sample that likely would have dissolved in the evolutionary time frame given. An evenly 
distributed binomial random variable with two possible outcomes (50/50) was used to simulate 

^1/2- 
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Lastly we chose a simple ratio of survival times to linearly scale the probability mechanism. 
Binaries with a separation greater than the 50% likelihood survival time will have a ratio greater 
than one and a better chance of dissociating (i.e. 50/50 turns into 75/25), while binaries with 
a separation less than the half point limit, will have a ratio less than one and a worse chance of 
dissociating (i.e. 50/50 turns into 25/75). For example: if a binary has a separation of W^A.U. 
its half-life is 3655Myr. If a binary evolves for 3655Myr, there will be a 50% chance of it being 
selected to dissolve. If the binary evolves for 10, OOOMyrs then: = 0.37 and 50 * 0.37 ~ 18. 
Therefore instead of a 50% chance of our binary surviving, there is an 18% chance of it surviving. 
Binaries that have been selected using our binomal random number generator, are removed in the 
posterior population samples. 

We note again that the dissolution algorithm is not a continuous function; its purpose is to 
generate in the final sample the expected number of dissolved binaries in a given population. The 
linear estimation of what should be an exponential (Poisson) function, will cause the dissolution 
algorithm to under-sample the number of dissolved binaries in our distribution. The dissolve or 
survive approach also does not allow for external forces and perturbations to expand the binaries. 
However, as stated above, the residual term from our model vs. observation comparison will 
provide some insight to the degree of which the sum of these external forces will effect our model. 
Binaries that have been removed are not considered in the following analysis. 



4 Results 

Rather then modeling individual binaries, we attempted to produce theoretical separation fre- 
quency distributions to compare with observational statistics. Histogram bin sizes were chosen to 
be similar for like quantities across the various cases. We then transformed our model parameters 



into model observations and compared our results with the observed Luyten (1969) separation dis- 
tribution. The statistical characteristics of the distribution of primary masses, secondary masses, 
separations and periods for the evolved binaries is given in Table [2] 
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4.1 Mass Graphs 

For comparison purposes, we plotted orbital frequency distributions for both the initial and pos- 
terior populations. Figure [T] and Figure [2] include not only the initial burst binaries but also those 
born in bursts later on. Thus, not all binaries formed at to = 10^ in the initial sample, while the 
posterior populations have identical t final- 

Figure [T] and Figure |2] are a representation of the same sample at an age of ~ lOGyr. In Figure 
|2| the random q-ratio has produced, in the secondaries, a distribution similar to the primaries. 
There are two artifacts in our secondary mass population that may not be physically realistic. 



First is the cutoff at O.OTM©, resulting from adopting the limit suggested by Chabrier (2003). 



However, two points mitigate this arbitrary cutoff: (1) as shown by Chabrier, the mechanism 
of star formation does have a minimum mass, and (2) imposing our q values to be a number 
between and 1 causes the average secondary mass value to be less than the primary, as expected. 



According to Chabrier the initial mass function at lower masses is poorly determined. Moreover, 
the LDS and SDSS surveys have detected relatively few M, L, and T dwarfs, particularly in binary 
systems, and not yet enough to warrant an attempt to model this region. 

A look at the posterior distributions in Figures [T] and [2] shows the impact of mass loss. Peaks 
in the distribution near the O.SMq to I.OMq range and a decrease of high mass stars is a sign that 
some of the MS stars have evolved to less massive WD stars, or Neutron Stars and Black Holes 
(NS and BH). At lOGyr, in some pairs both components have evolved to ~ O.GMq (WD + WD). 

Table 2: Statistical Characteristics of Evolved Distributions. 



Metric 


Average 


StdDev 


Skew 


Kurtosis 


Min 


Max 


Mass 1 (logMo) 


-0.36 


0.30 


-0.40 


-0.36 


-1.10 


1.01 


Mass 2 (log Mo) 


-0.58 


0.33 


0.27 


10.990 


-1.10 


0.81 


Separation log A.U. 


2.97 


0.55 


0.069 


-1.04 


1.74 


4.88 


Period log Myr 


-1.48 


0.84 


0.075 


-1.02 


-3.60 


1.11 
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4.2 Mass Ratio 

Figure [3] is a histogram of the mass ratio distribution for the initial population of binaries and 
the posterior population of binaries. There is a trend towards g ~ 1 in the initial population, as 
expected from the mass ratio generating function. A unique by-product of evolution is evident in 
the posterior population of the g-distribution: high mass MS stars have evolved into lower mass 
WDs. Now the secondary has a higher mass than the WD, resulting in a distribution of binaries 
with g > 1. This time, however, the g > 1 range represents the region where the primary has 
undergone PMS mass-loss and now has less mass than the secondary, i.e. one or both of the 
components of the binary have evolved into degenerate stars (WD, NS, or BH). The additional 
evolution to WD of the secondary would likely only cause a density change in the population 
above g = 1 or in the compact region near g = 0.9. 

4.3 Separation and Period 

Figures |4] and [5] present comparisons of the initial and posterior populations of binaries, the 
separations (in log {A.U.)) and the period distributions (in log (Myrs)). The sharp observed initial 
distribution of separations is a result of the initial parameters: a cutoff at lO'^A.U. (GMC birth 
cloud), and a lower cutoff (lOOA.U.) based on the definition of what qualifies as a FB. We impose 
an Opik distribution on top of this to produce the desired frequency distribution. The scatter 
is a result of binning and low number density. At higher densities (10^ binaries and greater), 
this scatter becomes negligible. Our initial distributions are leptokurtic (i.e. fiat-top), which is 
expected in log (a) due to our initial assumptions of a uniform log distribution. The skewness 
of the initial distribution ~ results from the flat sides and little or no discernable tail of the 
distribution. We see similar characteristics in the initial log (P) distributions(Figure |5|, however 
the edges of the distribution have been rounded out resulting from the effect the mass distribution 
has on the generated distribution. Our initial limits in the separation distribution results in the 
min (a > lO'^A.U.) and max (a < lO'^A.U.), respectively. 

Comparisons to posterior distributions in Figure |4] and Figure |5] demonstrate that the maximum 
separation is now log{A.U.) = 4.88 (75, 850^4. f/.). Two other parameters have also changed: the 
shift to higher separations has depleted the number of close binaries and increased the higher sep- 
aration binaries. This has effectively rounded out the histogram of separtions (kurtosis = —1.04) 
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and narrowed the histogram, (standard deviation = 0.55) compared to the standard deviation of 
0.58 in the initial distribution of Figure |4| Figure [5] shows a similar trend in variance, average and 
kurtosis. The period has the addition of a randomly generated mass factored into its calculation; 
this in turn has affected the kurtosis and skewness which now reflect a more normal distribution, 
per the central limiting theorem (CLT), as well as making the curve appear smoother. 

4.4 Analysis by Binary Category 

Figures |6] to [9] are overlay plots of the three binary classes MS + M S , WD + MS and WD + WD. 
All distrbutions are derivitives of the posterior population distribution described in the previous 
section; the "other" catagory refers to binaries with either component not in the MS or WD phase 
(RG, BH, NS...etc), the breakdown of population distribution types is given in Table |3] The 
statistical properties of each of the three samples are given in Tables |4] - [6j 

Table 3: Distribution of Evolved Binaries in Posterior Population. (Population size is 90,523). 



Metric 


MS + MS 


WD + MS 


WD + WD 


Other 


Count 


76,886 


7,327 


4,791 


1,519 


Percent 


84.94% 


8.09% 


5.29% 


1.68% 



4.4.1 MS+MS 

The MS + MS parameters are relatively unchanged and appear almost identical to the initial 
distribution plots in g-ratio, log (Myr) and \og{A.U.). However, the MS + MS mass plots do 
show missing upper mass stars. These, of course, have been transferred to the other groups as WD, 
BH, NS, or other stars. Similarly, a good number of higher separation binaries have dissolved, 
resulting in the observed change in skewness (observed asymmetry in separation distribution). Of 
the given model (90,523 binaries), 76,886 (84.94%) binaries remain in the MS + MS category. 
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Figure 6: Histogram of the distribution of primary masses from the MS+MS populaton (in red 
with 45 deg. hatching), from the WD+MS population (in green with 135 deg. hatching) and the 
WD+WD population (in blue with deg. hatching). Note the mixed colors represent overlaping 
distributions, bin sizes are in 0.05 log(M0) intervals. 

Table 4: Statistical Characteristics of Evolved MS+MS Distributions 



Metric 


Average 


StdDev 


Skew 


Kurtosis 


Min 


Max 


log a{A.U.) 

log P{Myr) 


2.94 
-1.49 


0.553 
0.841 


0.072 
0.066 


-1.12 
-1.05 


2.00 
-3.28 


4.00 
0.39 



4.4.2 WD+MS 

We observe the peak mass distribution of the primary in the WD + MS (i.e., the WD) at slightly 
less then unity {WD with masses of less then one solar mass). In addition, the secondary is 
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Figure 7: Histogram of the distribution of mass ratios from the MS+MS populaton (in red with 
45 deg. hatching), from the WD+MS population (in green with 135 deg. hatching) and the 
WD+WD population (in blue with deg. hatching). Note the mixed colors represent overlaping 
distributions, bin sizes are in 0.05 intervals. 

by definition an MS star (has not evolved yet). The secondary mass distribution has a peak at 
~ O.Tq, the approximate turn-off-mass for Disk stars, thus restricting the g-ratio distribution. The 
histogram for WD+MS g-ratio values (Figure]?]) shows a number of g > value binaries, reflecting 
the fact that some original secondaries are now more massive than their evolved primaries. 

The WD + MS histogram of separations shows an increase in excess kurtosis, i.e., the move- 
ment of low period and small separations towards higher values. This causes the more normal 
distribution curve that is seen {kurtosis =^ 0). Additionally, the average separation has increased 
to log (A.U.) = 3.1 ± 0.5, i.e., WD + MS pairs display an average separation a = 477A.U. greater 
than the original MS + MS distribution log{A.U.) = 2.9 ± 0.5. Thus, orbital amplification due 
to PMS mass-loss is shown to have increased the average separation of the population. We note. 
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Figure 8: Histogram of the distribution of separation from the MS+MS populaton (in red with 
45 deg. hatching), from the WD+MS population (in green with 135 deg. hatching) and the 
WD+WD population (in blue with deg. hatching). Note the mixed colors represent overlaping 
distributions, bin sizes are in 0.05 log{A.U.) intervals. 

Table 5: Statistical Characteristics of Evolved WD+MS Distributions 



Metric 


Average 


StdDev 


Skew 


Kurtosis 


Min 


Max 


log a{A.U.) 

log P{Myr) 


3.13 
-1.40 


0.489 
0.735 


0.262 
0.246 


-0.90 
-0.89 


2.24 
-2.78 


4.48 
0.72 



however, that this is not the order of magnitude difference that Greenstein ( 1986[ ) expected. Per- 
haps our models have underestimated the effects of environmental perturbations (GMC, chance 
encounters. Galactic tidal forces). There are 7327 WD + MS binaries in the final sample (i.e.. 
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Figure 9: Histogram of the distribution of period from theMS+MS populaton (in red with 45 deg. 
hatching), from the WD+MS population (in green with 135 deg. hatching) and the WD+WD pop- 
ulation (in blue with deg. hatching). Note the mixed colors represent overlaping distributions, 
bin sizes are in 0.1 log (Myr) intervals. 

~ 8.1%). 

4.4.3 WD+WD 



To date, there are only about two dozen (Sion et al. , 1991) double degenerate systems known 



and certainly not as many as shown here (A^ = 4791, ~ 5.3%). The most plausible reason for this 
lack of discovery of WD + WD binaries results from either one or both of the components being 
too dim to be detected in the surveys currently available. The histogram of the g-ratio (Figure [?]) 
show a tendency to nearly identical primary and secondary masses. This is to be expected with 
the narrow range of white dwarf final masses, resulting from the SSE's IFMR calculations. 

As predicted, the secondary companion (lower initial mass) produces a lower mass WD than 
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Table 6: Statistical Characteristics of Evolved WD+WD Distributions 



Metric 


Average 


StdDev 


Skew 


Kurtosis 


Min 


Max 


log a{A.U.) 
log P{Myr) 


3.15 
-1.36 


0.500 
0.749 


0.316 
0.308 


-0.81 
-0.82 


2.27 
-2.74 


4.66 
0.9 



the primary, causing the peak at g = 1 (Figure [Tj Similar to the WD + MS histogram of 
separations, the WD + WD histogram shows an increase in kurtosis, i.e. loss of short period and 
small separations, evolving towards higher values. This causes the more normal distribution curve 
that is seen {kurtosis =^ 0) compared to the MS-I-MS population. The increase is the result 
of continued evolution of the binaries. The post-MS sudden mass-loss of the binary secondary 
influences the separation and period distributions a second time, pushing the uniform distribution 
parameters closer to a more random, Gaussian distribution. Additionally, average separation has 
increased from \og{A.U.) = 3.13 for the WD + MS, to \og{A.U.) = 3.15 for the WD + WD, 
i.e. an average separation increase of 63.6A.U.. We also note that the maximum separation has 
increased to log{A.U.) = 4.66 range, much greater than the initial log{A.U.) = 4.0 maximum 
imposed on the initial distribution. 



4.5 Comparison to Observed Data 

We begin our comparison by testing a variety of distance models against the LDS WD + MS 
sample. For our models, we used a lO'^Myr evolutionary time with 10^ initial binaries (at t = 



0) and a ^°Syr'' ^^^^^ ^^^^ (i-^- constant average birth rate. Hurley et al. 2000A|. For each 



component, we have roughly determined absolute ugriz magnitudes using the Kurucz (1979) 



models. We found temperatures and gravities from the SSE output or from Chabrier (2003) 
evolutionary models for solar metallicity low-mass stars. Using the luminosities and separations, 
we applied a random distance model derived by assuming a uniform synthetic Solar Neighborhood 
population to approximate angular separation and apparent magnitude. This yielded a theoretical 
population histogram, restricted by the physical observation limits of the LDS. The distribution 
of the subset of modelled observed evolved binaries is given in Table [T) Note that, m comparison 
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to Table [3} there are fewer MS + MS observed than in the simulated population, an effect likely 
caused by the dimness of some low mass companions. Similarly, the percentage of WD + MS 
pairs in the observed distribution is higher than in total posterior population. The latter is likely 
caused by the expected relative brightness of the binary pair; as mass ratio tends towards one, a 
pair with a WD as a primary is likely to have a secondary with mass just slightly less than the 
progenitor of the WD, resulting in a brighter pair on average, than observed in the MS + MS 
distribution. 

Table 7: Distribution of Theoretically Observed Evolved Binarys in Postieror Population. Total 
observed population size is 29,6335 



Metric 


MS + MS 


WD + MS 


WD + WD 


Other 


Count 


23,507 


3,801 


1,442 


885 


Percent 


79.3% 


12.8% 


4.9% 


3.0% 


Fraction of Posterior 


30.6% 


51.9% 


30.1% 


58.3% 



4.6 Distance Model Comparison 

The comparison of our model and observed angular separation distributions is shown in Figures 



11 and 12 We expect the observational data, which is a much smaller sample (~ 6124 binaries) 
to have more random scatter than the modeled data (WD+MS population = 1670 binaries). 
Although the agreement is good dispite the scatter in the LDS sample, there are discrepancies 
(that we intend to address in a later paper). 

• The kurtosis of the model is is smaller than in observed sample. Perhaps wide binaries at 
smaller separations were under-observed by Luyten. Alternatively, our model's lower limit 
on separation could be unrealistic. 

• The scatter in the observed data is la greater than the modeled data. 

• The mean and skewness of both data sets is similar; this suggests that our basic model is 
sound, at least first order approximation. 
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5 Discussion 



The distribution of binary angular separations in our model is similar to the observed distribu- 
tion(i.e., the LDS sample). The LDS sample has a lognormal separation distribution that most 
likely stems from a log-uniform distance distribution. The estimated range of the LDS sample is 
between 1 < logd < 3 (where distance is in parsecs), similar to our model distance range. Initially 
we can compare the precentage of WD + MS pairs observed in the LDS sample (511/6124 = 0.08) 
with the precentage of modeled WD + MS pairs (.12) as well as the population of WD + WD pairs 
observed in the LDS sample (24/6124 = 0.003) with the population of WD + WD modeled (0.04). 
These distribution numbers suggest either we have underestimated the number of binaries lost to 
dissolution, or over-estimated the brightiness of some of the pairs. Likewise, these comparisons 
may also indicate where the LDS survey is under-sampled (i.e., the WD+WD pairs). 



Figures 10 and 12, present the results from our comparison of MS + MS LDS and model 
pairs as well as WD + MS LDS and model pairs. The figures show the cumulative distribution 
functions of the populations. Note that the separations are given in log(s") and have been put into 
0.2 interval bins. Our initial hypothesis that PMS mass-loss significantly distorts a distribution of 
separations is apparent in Figure [TT| the offset in the CDF curves reflects an increase in average 
angular separation. 

We use the Cramer- Von Mises goodness-of-fit test for comparison of the two-population cumu- 



lative distributions (Anderson, 1962), given in equations 18 and equation 19 Here, goodness-of-fit 
is defined as the hypothesis Hq : f = 0.0 where f = Fn{x) — Gm{,x), where F/v(x) and Gm{,x) are 
the cumulative distribution functions of samples of size and M, respectively. The Cramer- Von 
Mises goodness-of-fit statistic here is defined by T. 

NM r°° 

T = J [F^ix] - GM{x)fdH^^M{x) (18) 

NM ^ ^ 

For the MS + MS model and observation population set the values are = 26 and M = 26. 
For the WD + MS model and observation population set the values are also A^ = 26 and M = 26 
as we have used the same spacing for both empirical distribution functions. These empirical 



28 




Figure 10: Comparison of observed angular separation CDF for MS+MS models and WD+MS 
models. Red Dashed Line: MS+MS, Blue Solid Line: WD+MS 

distribution functions are based on the cumulative distributions generated for the populations 
with a bin interval of 0.2 between 0.4 < log (s") < 2.5. We compute T values for the comparison 
of MS+MS populations (T = 0.1343) and the comparison of WD+MS populations (T = 0.0418). 
Based on these statistics we estimate the limiting distribution, ai(T), of T for our models as in 



Anderson and Darling (1952): 



ai(T) = lim Pr{nuP < T} 



(20) 



For the two comparisons, from the Anderson and Darling (1952) tables aMs+MS > 0.5 and 



awD+MS < 0.5, i.e., the probability of the MS+MS model being drawn from the same distribution 
as the MS + MS observed population is greater than 10%, while the probability of the WD + MS 
model being drawn from the same distribution as the WD+MS observed population is at least 5%. 
Given our original hypothesis Hq : f = 0.0, the MS+MS distribution comparison hypothesis fails. 
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Figure 11: Plot of CDF of Model to Observed Cumulative Angular Separation Distributions for 
MS+MS pairs. Red Dashed Line: Observed, Blue Solid Line: Model 

and the WD + MS comparison hypothesis is accepted. We conclude then, that the distributions 
of observed and modeled WD + MS angular separations are similar, but not for the MS + MS 
distributions. 

We further demonstrate the relationship between the models by graphing the residual of the 



angular separation distributions in Figure [I3| The residuals plotted have no particular linear trend, 
nor offset. This suggests either a difference in spread or mean. There are however, offsets in the 
MS + MS comparison. These offsets cause the failure of the null hypothesis for the MS + MS 
comparison. This suggests: (1) there are additional factors for which we have not taken into 
consideration (such as inclination) and/or (2) small systematic errors associated with our model. 
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Figure 12: Plot of CDF of Model to Observed Cumulative Angular Separation Distributions for 
WD+MS pairs. Red Dashed Line: Observed, Blue Solid Line: Model 

6 Conclusions 

We have shown that PMS mass loss significantly shifts and reshapes the distribution of MS + MS 
separations. The semi-major axis of a wide binary is a difficult parameter to determine, requiring 
distance, inclination, period and other measurements. It is our hypothesis that a distribution of 
angular separations, preserves information about the effect of PMS main-losses on the binary 
system. We take as proof, the following observation: Both our MS + MS and WD + MS models 
compare reasonably well to the observed LDS distributions. Our model, however, has only been 
affected by PMS mass loss and a crude attempt has been made to account for other disrupting 
factors. This assumption and other first order approximation in the creation of the model, likely 
affect the goodness-of-fit test for the MS + MS observed populations. 



Figures 11 and 12 show the cumulative distribution of log (s") in the MS + MS and WD + MS 
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Figure 13: Plot of the Residual Square Model to Observed Cumulative Angular Separation Dis- 
tributions for MS+MS pairs (in blue open circles) and WD+MS pairs (in red closed circles) 

pairs. The difference in observational and theoretical distributions is apparent. This difference is 
a function of both the incompleteness of the survey and secondary effects for which we have not 
taken into consideration (such as inclination or external perturbations). For example, observations 
and models have shown that metallicity may play a much larger role in the amount of post Main- 



Sequence mass-loss (see Zhao et al. 2011B). 



Future work should include a more realistic model for stellar formation history rather than 
the constant birth function implemented here. Furthermore, this stellar birth model should be 
used to drive the initial conditions of both the metallicity and the mass distribution function. In 



view of the trends in the metallicity, especially as a function of stellar age (Zhao et al. , 2011B), 
additional effect on the mass-loss of the population as a function of time would be expected. 
We have demonstrated in this paper, and it has been demonstrated in other simulations, that 
post-MS mass-loss on average results in an increase in semi-major axis for visual and fragile 
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binaries. Changes to the metaUicity history would be expected to influence orbital amplification 
as a function of time, in turn affecting the final angular separation distribution. Our photometric 
colors were based on a consistent metaUicity for the population of binaries. Adjustments to the 
metaUicity history would affect the associated color indices, as well as the apparent magnitude of 
the binary components. Again, this would effect the visual angular separation distribution. 

The observed MS + MS and WD + MS distributions are affected by observational bias and/or 
physical orbital amplification. If orbital amplification results in an increase in orbitla separation, 
we conclude that binaries with small separations have expanded to middle separation binaries, 
middle binary separations have increased slightly, and binaries that were already wide may have 
been lost to dissociation. Our current simplistic model describes the orbital frequency distribution 
of FBs reasonably well. In general, our biggest issue is population size of the observed distribution. 
The statistically small observed population causes a large scatter in the frequency distribution, or 
rather a larger scatter than in the model distribution. 

The DS code is currently deficient at both extremes of stellar mass ( i.e., M < 0.07Mq and 
M > 9OM0 ). We do not have an evolutionary code applicable to these regimes, nor does our code 
include metaUicity variation. However, such variance is not of major important in the modeling 
of disk stars. Future work will require such an extension to the SSE code and extrapolation of 
the Kurucz color /atmosphere models. Low-mass models are available and have been taken into 



account for masses as low as M ~ O.OTSMq (Baraffe et al. , 1998). While their contributions to 
the current observed sample of FBs are small, in order to produce a realistic model of the evolved 
population, they should be taken into account. This lower mass limit does not directly affect the 
current project. 
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